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Abstract. - We study oscillation death (OD) in a well-known coupled-oscillator system that has 
been used to model cardiovascular phenomena. We derive exact analytic conditions that allow 
the prediction of OD through the two known bifurcation routes, in the same model, and for 
different numbers of coupled oscillators. Our exact analytic results enable us to generalize OD as 
a multiparameter-sensitive phenomenon. It can be induced, not only by changes in couplings, but 
also by changes in the oscillator frequencies or amplitudes. We observe synchronization transitions 
as a function of coupling and confirm the robustness of the phenomena in the presence of noise. 
Numerical and analogue simulations are in good agreement with the theory. 



Coupled oscillator systems exhibit a variety of phenom- 
ena relevant to physics, biology, and other branches of 
science and technology. Here, we study oscillation death 
(OD) [1], a form of synchronization [2] in which the oscil- 
lators interact in such a way as to quench each other's os- 
cillations [3-6] . This intriguing phenomenon was noted in 
the 19th century by Rayleigh [2], who found that adjacent 
organ pipes of the same pitch can reduce each other to si- 
lence. Since then, OD has been studied in diverse applica- 
tions including oceanography [7] , chemical engineering [8] , 
solid-state lasers [9] and a variety of other experimental 
systems [4,10-12]. OD is known to occur via two distinct 
bifurcation mechanisms: (i) Hopf bifurcation, where the 
coupling induces stability at the origin of the phase space, 
thus collapsing the orbits to zero, which can happen only 
if the oscillators are sufficiently different [5, 13-16] (or for 
identical oscillators if there are delays [17, 18] in the cou- 
pling); or (ii) for non-identical oscillators, saddle- node bi- 
furcation [4] in which new fixed points appear on/near the 
coupled limit cycles, annihilating the periodic orbits. Re- 
cently, Karnatak et al. [19] were able to produce OD in 
two identical coupled oscillators through the saddle-node 
route, using dissimilar non-delayed coupling. 

In this Letter, we show that a coupled-oscillator sys- 
tem, which has been used extensively in modeling cou- 
pled rhythmic processes in mathematics, physics and biol- 
ogy [27] can undergo OD via both bifurcation routes. We 
obtain exact analytic conditions for OD, and compare the 
theory with numerical simulations and analogue electronic 



experiments. We thus generalize OD as a phenomenon 
that occurs, not only through coupling-increased dissipa- 
tivity [2] , but also when a measure of dispersion among the 
parameters of the coupled system is exceeded [21]. We also 
show that, near the onset of death, the coupled system al- 
ternates between periodic, quasiperiodic and even chaotic 
behavior, reflecting the complex temporal variability ob- 
served in real biological systems [22,23]. 

Our model is a set of five coupled oscillators that suc- 
cessfully reproduces many phenomena seen in the cardio- 
vascular system (CVS), e.g. modulation [20] and synchro- 
nization [24]. Each oscillator has its own characteristic 
frequency and amplitude [25] (see Table 1) and emulates 
a particular physiological function. Heart and respiration 
are obvious physiological processes; the myogenic oscilla- 
tion is related to the intrinsic self-regulatory activity of the 
smooth muscle tissue in the walls of the blood vessels; the 
neurogenic oscillation is associated with the neural control 
by the central nervous system; and the nitric oxide related 
endothelial oscillation is associated with metabolic activ- 
ity mediated by the endothelial tissue that lines the whole 
CVS [26]. 

The basic unit. 



Xi = -XiQi - yiCJi + P^(X, Y) 
Vi = -ViQi^ XiUi + (X, Y) , 



(1) 



is the Poincare oscillator [27] with: i = 1, . . . , m, (m = 5); 



211 fi and a are con- 



Suarez- Vargas, Gonzalez, Stefanovska, McClintock 



Table 1: Typical values of frequency (Hz) and relative ampli- 
tude (arbitrary units) of physiological rhythms in humans, as 
measured in blood flow by laser-Doppler flowmetry and anal- 
ysed by wavelet transform. The characteristic frequencies vary 
around the values indicated; the amplitudes are only estimated 
based on measurements [20,25]. 



Activity 

Heart 

Respiration 
Myogenic 
Neurogenic 
Endothelial 



Frequency 
1.1 Hz 
0.3 Hz 
0.1 Hz 
0.04 Hz 
0.01 Hz 



Amplitude 
0.5 a.u. 
1.0 a.u. 
1.0 a.u. 
1.0 a.u. 
0.5 a.u. 




stants that represent amplitude, frequency and dissipation 
rate respectively. X and Y G ^ and Pi^Qi : R 
are scalar coupling functions. We suppose that the inter- 
oscillator interactions can be approximated by a mean 
field, so that Pi(X, Y) = ^ ^^j, and Qi(X, Y) = 0. 

We now review briefly some basic concepts, assuming 
an autonomous n-dimensional dynamical system: 



F(Z), 



(2) 



where Z G R^ and F : R^ R^ is a general nonlinear 
vector function. All the points, Z*, in phase space satis- 
fying the equation F(Z*) = 0, are called fixed points of 
the dynamical system (2). Their stability and quality can 
be determined by the eigenvalues of the Jacobi matrix, 
provided their real parts are nonzero: 



det 



' aF(Z ), 
dZ 







(3) 



We consider first the case of two coupled oscillators, mod- 
elling e.g. the cardio-respiratory interactions. 



Xi 

m 



-Xi 



a{^/xl 



= -yia{^xlTyl 



-X20l{^J~^- 



ai) + x\ijj\, 

-0.2) -^2^2 + e(^i +^2), 

+ X2UJ2' 



(4) 

Here the origin is always a fixed point. For e = 0, this 
point is unstable (focus or node). Besides this fixed point, 
there is a stable limit cycle, corresponding to autonomous 
oscillations. Eq. (3) for the point Z* = and e 7^ in the 
dynamical system (4) will be: 



{aai + e) 

e 




A 



aai — A 





{aa2 + e) — A 



002 






—002 
aa2 — A 



= 0.1 



(5) 



From (5) we determine that the origin can become a stable 
fixed point if 

2aai + e < 0, 
2aa2 + e < 0, 



(6) 



Fig. 1: (Color online) Saddle- node bifurcation diagram showing 
the appearance of new equilibria and their stability for coupled 
system (4) . The dashed (red) line corresponds to unstable fixed 
points, while the continuous (blue) line depicts stable ones. 



in which case the limit cycle no longer exists. Thus, if 
a2 > ai, then e < — 2aa2 is a sufficient condition for OD. 
This is the Hopf bifurcation route to OD, and it occurs 
only for negative e. 

For e > 0, there is a critical value Cc such that, for 
e > ec, four new fixed points appear: two unstable and 
two asymptotically stable. In order to make numerical 
estimations from our model we make use of the extensive 
earlier research on the CV system [20, 25] embodied in 
the frequencies and amplitude relationships summarised 
in Table 1, so that we take ai/a^ ~ 1, ai/a^ ~ 0.5, for 
i = {2,3,4}. 

Fig. 1 shows the bifurcation diagram for xi , calculated 
for ai = 0.5, a2 = 1, /i = 1.1, /2 = 0.3 and a = 1. Note 
that for this range of e > the origin is always an unstable 
fixed point. The end of the oscillations is marked by the 
appearance of the new fixed points at e ~ 3.578. They are 
obtained from the set of algebraic equations: 



-xia{^ yxl + yf - ai) - yiui + e{xi + X2) 
-yxOL{^Jx\ ^y\ - ax) + x^ijJx 



0, 
0, 



-X2Ci{^ Jx\ ^y\ - a2) - ^2^2 + e{xi + X2) 
-?/2a( V^fT^ - a2) + X2UJ2 

After some algebra, we obtain the relations: 

uji{xl^yl) = eyi{xi^X2), 
^2(^2 + ^2) = (^y2{xi^X2). 



0. 



(7) 



(8) 



Analysis of Eq. (8) combined with Eq. (7) shows that, for 
e <C 2uJi and e <C 2cj2, the only fixed point is (0,0,0,0). 
For sufficiently large values of e, we obtain four new fixed 
points. The fixed points can appear only in pairs of stable- 
unstable points. This is the saddle-node bifurcation route 
to OD, which occurs only for positive e. 
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We now define xi = ri cos yi = ri sin ($i), X2 



r2Cos(<l>2), y2 = r2sin(<l>2), where rf = xf 



-y2- As in our model cji ^ C(;2, <^2 > which implies 



[xi) <C (^2) • From Eq. (8) we obtain the equality: 
002 = ecos (<l>2) sin {^2)' 



(9) 



This is equivalent to the equation: 



2002 

e 



sin(2<l>2). 



(10) 



As sin(2<l>2) < 1, this means that we have new solutions 
to the algebraic equations (7) only when 



e > 2uo2. 



(11) 



Eq. (11) provides a simple and understandable analytic 
estimate of the critical value for the bifurcation. 

Fig. 2 (top) illustrates a numerical simulation of the 
coupled oscillators (4) near criticality, showing how the 
oscillations evolve as e is increased by steps of 0.1. We 
have established that, for the chosen parameter values, 
the oscillations die when Cc ~ 3.578, gratifyingly close 
to the ec ^ 47r/2 = 3.77 predicted by (11). The differ- 
ence is attributable to the approximations made in de- 
riving Eq. (11). Using the same values of parameters, we 
also performed numerical simulations for negative e. Fig. 2 
(bottom), and found that OD occurs through supercritical 
Hopf bifurcation when e ~ —1.88, in agreement with the 
theoretical prediction (6). 

We have also used the above mathematical tools to in- 
vestigate the full model [20] with five coupled oscillators 
(z = 1, . . . , 5). Solutions of the algebraic equations 



= -XiQi - yji^i + e^^Xj 
= -yiqi^XiUJi, 



(12) 



correspond to the fixed points of (1) for five coupled oscil- 
lators. After some calculations we get the equations 



(13) 



We have found that, for e > ec, where Cc is some critical 
value, the system possesses four additional fixed points 
(two stable and two unstable). 

Because coi = 27r/^, and using the empirical physiologi- 
cal values (table 1) for the /^'s and a^'s mentioned above, 
we can use the inequalities uo^ < uo^ < uos < 002 < ^1^ 



and a( = 



in order to obtain 



an analytic expression for ec- The resulting approxi- 
mate equations lead to the following relations for the sys- 
tems' fixed points: XI ^ ujiai; x^ ^ y^ ^ y^ ^ 0; 

Xr^ ^ X^ ^ X^ X^ , where x^ 

yi ^ ai; t/^ = ^^^^^ ^ ^ 

The new simplified equation (13), with i = 2, is 



2^ ' 



(14) 



w 2 




r -1 




1000 1500 
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Fig. 2: (Color online) Top: Time series of xi(t) (bold curve) 
and X2{t) from numerical simulations of the system (4), show- 
ing saddle- node OD. Parameter values were the same as in 
Fig. 1. The dashed line shows the step-wise variation of the 
coupling constant e in a positive range. Note that after OD has 
occurred lim {xi (t), X2(t)} > 0. Bottom: Time series of xi{t) 

t — >oo 

and X2{t) from numerical simulations of the system (4), show- 
ing supercritical Hopf bifurcation OD. Parameter values were 
the same as in Fig. 1. The diagonal stepped line starting at the 
origin indicates how e was varied in the negative range. In this 
case lim {xi (t), X2(t)} = 0, after OD has occurred. Note the 

difference in abscissa timescales between the top and bottom 
parts of the figure. 



The new stable fixed points are possible when 
y^e^ + I2x^u02e - Aco^y^ > 0. 



(15) 



The asymptotically stable equilibrium points are the new 
attractors of the dynamical system. 

Using physiologically estimated values for the parame- 
ters (table 1) in (15) we obtain that Cc ~ 0.468. This result 
is in agreement with the numerical simulations shown on 
Fig. 3(a), where there is OD for e ^ 0.47. 

We have also modeled the dynamics of the five-oscillator 
system for the given parameter values using analogue elec- 
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time(s) 

Fig. 3: (Color online) Simulations of the coupled system (1) us- 
ing five oscillators, as e is increased in steps of 0.1. Parameter 
values are based on physiology (see text), (a) Numerical simu- 
lation showing the outputs from all five oscillators, exhibiting 
OD for e ^ 0.47. (b) Analogue electronic simulation, showing 
the output from the first oscillator xi. 




Fig. 4: (Color online) Behavior of the coupled system (4) for 
e — 3.578, very close to death. Row (a): Attractor in xi-yi 
space, time series xi(t) and its Fourier transform. Row (b): 
First-return maps of the instantaneous period and amplitude 
signals. 

tronic circuits [28]. The result of Fig. 3(b) shows the ex- 
perimental observation of OD in good qualitative agree- 
ment with both the theory and numerical simulations. 

We have found that, just before the onset of death, the 
2-oscillator system (4) exhibits highly complex, quasiperi- 
odic and chaotic, behavior. Fig. 4(a) plots the xi-yi phase 
diagram when e = 3.578 ~ Cc. For comparison the circle in 
the background shows the limit cycle when e = 0, and the 
two full (green) spots signal the place where the new fixed 



points will arise. The attractor's structure, the seemingly 
random amplitudes of the time-series, the lengthening of 
the period with a square-root scaling, and the power-law- 
like spectrum, indicate that the system enters a chaotic 
regime just before it dies. Fig. 4(b) shows the first-return 
maps (FRM) of the discrete instantaneous period and am- 
plitude signals obtained by sampling xi{t) at the peaks. 
We observe strong period and amplitude variability in the 
FRM, with the period showing a multimodal distribution. 
The amplitude of peaks has a more uniform distribution 
showing fewer preferred values, and the amplitude differ- 
ence between successive peaks is distributed normally. In 
terms of the cardiovascular analogue this complex behav- 
ior might be seen as a predictor of imminent death of the 
cardio-respiratory coupled oscillations. Moreover, multi- 
modal period distributions are typical of the time variabil- 
ity observed in other biological systems, e.g. the intermi- 
totic time of human skeletal cells, and moments of change 
in the rotation direction of flagella [23]. It is evident that, 
at least in this kind of coupled biological system, the vari- 
ability of the oscillations increases near the onset of global 
bifurcations. 

Thus, OD in (1) can occur via both of the known routes. 
When e is negative and below a critical value, the origin 
becomes asymptotically stable and the oscillations die at 
almost constant frequency, i.e. the Hopf bifurcation sce- 
nario. In the second mechanism, OD occurs when new 
fixed points appear on the former attractor after e has sur- 
passed a positive threshold: the amplitude remains almost 
constant while the frequency decreases with a square-root 
scaling in a saddle-node bifurcation. 

When a system arrives in the OD regime it lies quies- 
cent. Although apparently trivial, this state results from 
diverse complex interactions between the coupled elements 
that form the system, as we show in this Letter. Thus OD 
can be seen as another kind of complex collective motion, 
much in the same way as synchronization arises in complex 
coupled systems as a self-organized dynamics. 

In order to complete the analysis of system (1) in this 
context we investigate the relationship between synchro- 
nization and OD by measuring the synchronization in- 
dex [2], defined as: 71,1 = (cos ^1,1)^ + (sin ^1^1)^, where 
^1^1 is the relative phase difference between the oscilla- 
tors. This index measures quantitatively the strength of 
1:1 synchronization. 

Fig. 5 shows the calculated values of 71,1 as a function 
of e for the coupled system (1). Initially 71,1 is very small 
for low values of e and it increases monotonically in a 
quasi-exponential fashion as e is increased. However when 
e reaches a value near 3.15, 71,1 no longer grows monoton- 
ically, but has alternating epochs of growth and decrease, 
corresponding to topological changes in the structure of 
the attractors of both oscillators that lead to complex 
transitions in the 1:1 synchronization state. These results 
suggest that the transient epochs of cardio-respiratory 
synchronization seen (for higher synchronization ratios) in 
many studies of resting humans, both awake [29-31] and 
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Fig. 5: Synchronization index 71,1 as a function of the coupling 
strength e for the system (1). 

asleep [32], may arise in part from changes in coupling, 
as well as from drifts in the natural frequencies and other 
parameters. 

Finally, in order to explore the robustness of our results 
in a more realistic framework, we performed simulations of 
system (1) including the presence of stochastic forces, since 
they can lead to spurious detection of complex phenomena 
[33]. Our stochastic model is: 

Xi = -XiQi - yiUJi ^ e{xi ^ X2) ^ D^:jc 
Vi = -yiQi^XiiJi, 

where corresponds to white Gaussian noise with zero 
mean and variance D^. 

The simulation of Eq. (16) is shown in Fig. 6. The 
parameters used in this noisy equation are the same as 
those used for Eq. (1), besides D = 1.0. In the top panel, e 
was piecewise increased through a range of positive values, 
and OD was observed to occur at a value very close to the 
noiseless case, e = 3.57. Similarly in Fig. 6 (bottom) e 
was piecewise increased but now in the negative direction. 
OD was still observed, at a value near that of the noiseless 
case, e = —1.8. Thus the noise term in Eq. (16) does 
not eliminate, or change the nature of, the bifurcations 
leading to the appearance of OD; however, it influences the 
asymptotic value of the fixed points to which the system 
evolves. 

In our study of OD we have obtained analytic relations 
predicting the onset of OD in the CV model, both for two 
and five coupled oscillators, using assumptions provided 
by experimental physiological measurements. In addition, 
all the theoretical results have been verified in numerical 
simulations. We have also observed partial OD [34] in our 
numerical model, where some oscillators die while others 
remain active. It occurs if the eigenvalues corresponding 
to certain variables Xi and yi have negative real parts, 
while those corresponding to other oscillator variables do 
not have this property. For instance, this happens if 2aai-\- 



4 




Fig. 6: Temporal behavior of noisy system (16) as e is varied 
piecewisely. Top: e increases positively and OD is observed 
near e=3.57. Bottom: e decreases negatively and OD is ob- 
served near e=-1.8. In both cases noise intensity is D = 1. 

e < 0. See Fig. 2 (bottom). 

Discussions of OD usually focus on the strength of the 
coupling coefficient. The common phrase is that "for large 
couplings, amplitude death will take place" . However, now 
that we have analytic expressions for the critical value, 
also depending on other parameters, we can predict the 
onset of OD when different parameters are changed: OD 
is evidently a multi-parameter-sensitive phenomenon. It 
can be induced, not only by changes in couplings, but also 
by changes in the oscillator frequencies and amplitudes. 
For instance, the analytic conditions for OD given by (6) 
indicate that, if the amplitude parameters ai and a2 de- 
crease (in fact it is sufficient to decrease the largest one), 
then OD can occur for two coupled oscillators even when 
the coupling coefficient remains constant. The same can 
happen if the excitability parameter a is decreased be- 
low a critical value. Another interesting feature of OD 
inferred from (11) and (15) is that, if e > 0, and e > ec, 
new asymptotically stable fixed points can appear. Thus 
the oscillations die, but the dynamical variables can take 
non-zero asymptotic values. We remark that the critical 
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value Cc depends on the oscillators' frequencies. That is, 
even for fixed coupling strength e, if the frequency of one 
of the oscillators falls below some critical value, then all 
the coupled system can stop dead. This behavior is true 
both for two or five coupled oscillators. 

What are the implications of these results for the CVS? 
Given that (1) successfully models many features and 
states of the CVS, we may speculate that phenomena anal- 
ogous to OD may also occur there. It would mean that, 
for certain parameter combinations, the mutual interac- 
tions between the cardiac, respiratory, myogenic, neuro- 
genic and endothelial oscillators might serve to bring one 
or all of them to a halt - or at least to inhibit function. 
Treatment to prevent or remedy such a scenario might in- 
volve the use of drugs to modify some of the parameters a, 
ai or Ui so as to take the individual away from the regime 
of danger. 

These results could also be applicable to coupled ar- 
rays of neurons in the brain. In the first demonstration 
of OD in a biological system Ozden et. al [10] coupled a 
man-made device to an array of real neurons and OD was 
provoked by taking the system into the strong-coupling 
regime. We hypothesize that similar results could be ob- 
tained in response to frequency changes, according to rules 
similar to (6) or (11). The system would not then need 
to be in the strong-coupling regime, which in some cases 
could damage sensitive biological tissue. 

In summary, we have shown analytically and through 
simulations that OD in the CVS model (1) can occur via 
either of the known bifurcation mechanisms. Furthermore, 
OD can be induced, not only by an increase in coupling 
strength as conventionally accepted, but also by changes 
of e.g. amplitude and/or frequency. We have shown it to 
be sensitive to many different parameters. Connections 
between OD in the model, and phenomena in the CVS, 
appear possible but remain to be explored. 
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